function [p, grid, Theta, Delta, Gamma, S, i] = PDEFD(sigma, r, K, S0, dx, dt, T, ...
   scheme, type, XminFactor, XmaxFactor, boundaryType)
%
% This function computes among others the premium of an option given a certain
% finite difference scheme.
% 
% 
% Arguments:
% PDEFD(sigma, r, k, s0, dx, dt, t, scheme, type, 
%     XminFactor, XmaxFactor, boundarytype)
% 
%  sigma       The volatility of the underlying stock.
%  r           The interest rate.
%  K           The strike price of the option.
%  S0          The price of the underlying stock at time 0.
%  dx          The horizontal step size.
%  dt          The vertical step size.
%  T           The maturity time of the option.
%  scheme      The finite difference scheme. Can be:
%                 FCTS  for forward centered time scheme.
%                 CNS   Crank-Nicolson scheme. 
%  type        The type of the option. Can be: put or call.
%  XminFactor  The minimum value of the horizontal axis.
%  XmaxFactor  The maximum value of the horizontal axis.
%  boundaryType   The type of boundaries XminFactor, XmaxFactor. Can be:
%                    Fixed    for boundaries which are not depedent of the
%                             horizontal step size
%                    Variable for boundaries which are dependen of the 
%                             horizontal step size.
%
% Output:
%  p           The approximation of the premium of the option.
%  grid        The grid which contains all the values which were used to
%              compute the premium. Note that it is correct for time i.e.
%              grid(:, 1) gives all the premiums at time 0.
%  Theta         The Theta of the option.
%  Delta       The Delta of the option.
%  Gamma       The Gamma of the optoin.
%  S           The values which were used on the horizontal axis.
%  i           The index in S and vt which gives the correct premium.
%

mu = (r - 0.5 * sigma^2);

if strcmp(scheme, 'FCTS')
   % Determine if the parameters are chosen such that the solution of the FCTS
   % will be stable.
   if (dt / dx^2 > 1/2 || dt / dx^2 <= 0)
      ME = MException('FCTS:instable', ...
         'please choose dx and dt such that 1 < dt / dx^2 <= 1 / 2');
      throw(ME);
   end
   % Set up the coefficients for the banded matrix which will be used in the 
   % iteration.
   a = [0; 2 * dx^2; 0];
   % Set up coefficients for the vector which will be used in the iteration.
   b = [mu * dt * dx + sigma^2 * dt; ...
      2 * (dx^2 - sigma^2 * dt - r * dt * dx^2); ...
      sigma^2 * dt - mu * dt * dx];
elseif strcmp(scheme, 'CNS')
   % Set up the coefficients for the banded matrix which will be used in the 
   % iteration.
   a = [-mu * dx * dt - sigma^2 * dt; ...
      4 * (1 + 0.5 * r * dt) * dx^2 + 2 * sigma^2 * dt; ...
      mu * dx * dt - sigma^2 * dt];
   % Set up coefficients for the vector which will be used in the iteration.
   b = [mu * dx * dt + sigma^2 * dt; ...
     4 * (1 - 0.5 * r * dt) * dx^2 - 2 * sigma^2 * dt; ...
     sigma^2 * dt - mu * dt * dx];
else
   % Notify the user that only the forward-time-centered - and Crank-Nicholson
   % scheme are implemented.
   ME = MException('PDEFD:wrongInput', 'only FCTS and CNS are implemented');
   throw(ME);
end

if mod(T / dt, 1) ~= 0
   ME = MException('PDEFD:wrongdt', ...
      'please choose dt such that T / dt is in N')
   throw(ME);
end

% Transform the initial price to the log coordinates.
X0 = log(S0);

% Use a function pointer for the payoff.
if strcmp(type, 'call')
   boundary = @(X, tb)BoundaryCall(X, K, r, T, tb);
elseif strcmp(type, 'put')
   boundary = @(X, tb)BoundaryPut(X, K, r, T, tb);
elseif strcmp(type, 'digital')
   boundary = @(X, tb)BoundaryDigital(X, K, r, T, tb);
else
   ME = MException('PDEFD:wrongInput', 'only call and put are implemented');
   throw(ME);
end

if strcmp(boundaryType, 'variable')
   % Set up the boundaries of the grid which will be used for the FD scheme.
   Xmin = X0 - dx * XminFactor;
   Xmax = X0 + dx * XmaxFactor;
elseif strcmp(boundaryType, 'fixed')
   Xmin = XminFactor;
   Xmax = XmaxFactor;
else
   ME = MException('PDEFD:wrongInput', ...
      'only variable and fixed boundary types are implemented');
   throw(ME);
end

% Determine the number of elements in the grid.
N = round((Xmax - Xmin) / dx);
% Determine the number of iterations.
M = T / dt;

% Initialize the diagonal/banded matrix.
A = diag(ones(N, 1) * a(2), 0)+ diag(ones(N - 1, 1) * a(3), -1) + ...
   diag(ones(N - 1, 1) * a(1), 1);

% Set up the initial vector which will hold the final pay off of the option
% given a number of stock prices.
v = zeros(N, 1);
S = zeros(N + 2, 1);
for i = 1 : N
   S(i + 1) = exp(Xmin + (Xmax - Xmin) / N * i);
   v(i) = boundary(exp(Xmin + (Xmax - Xmin) / N * i), 0);
end
S(1) = exp(Xmin);
S(N + 2) = exp(Xmax + dx);
% Initialize some arrays.
c = zeros(N, 1);
[n, m] = size(v);
grid = zeros(N + 2, M);

% Do a backwards iteration since tau walks back in time.
for tau = 1 : M
   % Set the boundary conditions.
   v0 = boundary(exp(Xmin), dt * tau);
   vN1 = boundary(exp(Xmax + dx), dt * tau);
   % Asstaugn the columns to the grtaud whtauch wtaull be returned.
   grid(2 : N + 1, M - tau + 2) = v(1 : N);
   grid(1, M - tau + 2) = v0;
   grid(N + 2, M - tau + 2) = vN1;
   % Set up the elements of the vectors.
   c(1) = b(1) * v(2) + b(2) * v(1) + b(3) * v0 - a(3) * v0;
   c(2 : N - 1) = b(1) * v(3 : N) + b(2) * v(2 : N - 1) + b(3) * v(1 : N - 2);
   c(N) = b(1) * vN1 + b(2) * v(N) + b(3) * v(N - 1) - a(1) * vN1;
   % Solve the systems of linear equations represented with the matrix A and
   % the vector c.
   v = A \ c;
end
% Set up last column of the grid which will be returned.
grid(2 : N + 1, 1) = v(1 : N);
grid(1, 1) = boundary(exp(Xmin), dt * (M + 1));
grid(N + 2, 1) = boundary(exp(Xmax + dx), dt * (M + 1));

for i = 1 : length(S)
   if (abs(S(i) - S0) < 0.00001)
      p = v(i - 1);
      break;
   end
end


Theta = zeros(length(grid(:, 1)), 1);
Theta(:) = (grid(:, 2) - grid(:, 1)) / dt;
Delta = zeros(length(grid(:, 1)) - 2, 1);
Gamma = zeros(length(grid(:, 1)) - 2, 1);

for j = 2 : length(grid(:, 1)) - 1
   Delta(j) = (grid(j + 1, 1) - grid(j - 1, 1)) / (S(j + 1) - S(j - 1));
   Gamma(j) = (grid(j + 1, 1) - 2 * grid(j, 1) + grid(j - 1, 1)) / ...
      (S(j + 1) - S(j - 1))^2;
end

% vim: ft=matlab:tw=79:et:autoindent
